Figure3S
Aina Rill-Hinarejos
2023-07-04
Introduction
This is the code to reproduce Figure 3S Additional longitudinal and repertoire analysis of T cells across timepoints.: Fig3S C-H. To obtain the data object used in this notebook, please run 01_TCR_Data_Analysis.Rmd.
Package loading
library(gdata)
library(Seurat)
library(reshape2)
library(ggplot2)
library(ggpubr)
library(cowplot)
library(scRepertoire)
library(circlize)
library(Startrac)
library(stringr)
library(plyr)
library(dplyr)
library(ggh4x)
library(ggbeeswarm)
library(monocle)Loading colors
cols_patient =c("Patient 1"= "aquamarine", "Patient 2"="lightpink", "Patient 3"="yellow1", "Patient 4"="skyblue", "Patient 5"="sienna1")
cols_patient2 <- c("aquamarine", "#6FDFBA", "lightpink", "#DF9FA9", "yellow1", "#DFDC00", "skyblue", "#76B4CF","sienna1", "#DF733E")
cols_timepoint<- c("IP"="#4E6AAB","Peak"="#e78ac3")
cols_anno <- c("CD4+ Naive T cells"= "#33A02C",
"CD4+ CEntral/Effector memory T cells (CM/EM)"="#B2DF8A",
"CD8+ cytotoxic T cells"="#185B88",
"CD8+ Effector T cells (E)"="#1F78B4",
"CD8+ Eff/Mem T cells (EM)"="#A6CEE3",
"Early prolif: MCM3/5/7+ PCNA+ T cells"="#FB9A99",
"Late prolif: histones enriched MKI67+ T cells"="#E31A1C",
"Late prolif: CCNB1/2+ CDK1+ T cells"="#CAB2D6",
"Late prolif: STMN1+ BIRC5+"="#FDBF6F",
"Ribosomal/Mitochondrial/Degraded cells"="#FF7F00",
"gamma-delta T cells"="#6A3D9A")
strip <- strip_themed(background_x = elem_list_rect(fill = c("#A6CEE3", "#1F78B4","#185B88", "#B2DF8A","#33A02C", "#CAB2D6", "#FB9A99", "#E31A1C", "#FDBF6F", "#FF7F00" )))
col_clono <- c("Hyperexpanded (100 < X <= 500)"="#810F7C", "Large (20 < X <= 100)"="#8856A7","Medium (5 < X <= 20)"= "#8C96C6","Small (1 < X <= 5)"= "#9EBCDA", "Single (0 < X <= 1)"="#BFD3E6", "No clonotype detected"="gray69")Data import
combined_infusion <- readRDS("../data/combined_files/combined_infusion_CARPOS.rds")
combined_infusion<- addVariable(combined_infusion, variable.name = "patient",
variables = c("patient1", "patient2","patient3", "patient4","patient5"))
table_combined_infusion<- clonalDiversity(combined_infusion,
cloneCall = "gene",
group.by = "sample",
n.boots = 100, exportTable = T) #diversity_infusioncombined_peak <- readRDS("../data/combined_files/combined_peak_CARPOS.rds")
combined_peak<- addVariable(combined_peak, variable.name = "patient",
variables = c("patient1", "patient2","patient3", "patient4","patient5"))
table_combined_peak<- clonalDiversity(combined_peak,
cloneCall = "gene",
group.by = "sample",
n.boots = 100, exportTable = T) #diversity_peaktable_combined_infusion$timepoint <- "Infusion"
table_combined_peak$timepoint <- "Peak"
table_combined_infusion$patient <- c("Patient 1", "Patient 2", "Patient 3", "Patient 4", "Patient 5")
table_combined_peak$patient <- c("Patient 1", "Patient 2", "Patient 3", "Patient 4", "Patient 5")
table_combined <- rbind(table_combined_infusion,table_combined_peak)
table_combined$sample <- factor(table_combined$sample)seurat <- readRDS("../data/MENENDEZ_DEF.rds")
seurat_carneg <- subset(seurat, subset = Class1 == "CAR-")
seurat_carpos <- subset(seurat, subset = Class1 == "CAR+")
seurat_carneg <- subset(seurat_carneg, subset = cloneSize != "No clonotype detected")
seurat_carneg <- subset(seurat_carneg, subset = cloneSize != "Single (0 < X <= 1)")
seurat_carneg_nogd <- subset(seurat_carneg, subset = annotation != "gamma-delta T cells")
seurat_carpos <- subset(seurat_carpos, subset = cloneSize != "No clonotype detected")
seurat_carpos <- subset(seurat_carpos, subset = cloneSize != "Single (0 < X <= 1)")
seurat_carpos_nogd <- subset(seurat_carpos, subset = annotation != "gamma-delta T cells")
seurat_carneg_IP <- subset(seurat_carneg_nogd, subset = Timepoint == "IP")
seurat_carneg_PEAK <- subset(seurat_carneg_nogd, subset = Timepoint == "Peak")data_pt1 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient1")
data_pt1_I <- subset(data_pt1, subset = Timepoint == "IP")
data_pt1_P <- subset(data_pt1, subset = Timepoint == "Peak")
data_pt2 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient2")
data_pt2_I <- subset(data_pt2, subset = Timepoint == "IP")
data_pt2_P <- subset(data_pt2, subset = Timepoint == "Peak")
data_pt3 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient3")
data_pt3_I <- subset(data_pt3, subset = Timepoint == "IP")
data_pt3_P <- subset(data_pt3, subset = Timepoint == "Peak")
data_pt4 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient4")
data_pt4_I <- subset(data_pt4, subset = Timepoint == "IP")
data_pt4_P <- subset(data_pt4, subset = Timepoint == "Peak")
data_pt5 <- subset(seurat_carpos_nogd, subset = Patient_id == "patient5")
data_pt5_I <- subset(data_pt5, subset = Timepoint == "IP")
data_pt5_P <- subset(data_pt5, subset = Timepoint == "Peak")data2<- readRDS("../data/dataordered_V2.rds")Plot Fig3S C
occRep_ip<-clonalOccupy(seurat_carneg_IP, x.axis = "annotation", exportTable =T)
occRep_peak <-clonalOccupy(seurat_carneg_PEAK, x.axis = "annotation", exportTable =T)
occRep_ip$Timepoint <- "IP"
occRep_peak$Timepoint <- "PEAK"
occRep <- rbind(occRep_ip, occRep_peak)
occRep$annotation <- factor(occRep$annotation, levels =
c("CD8+ Eff/Mem T cells (EM)",
"CD8+ Effector T cells (E)",
"CD8+ cytotoxic T cells",
"CD4+ CEntral/Effector memory T cells (CM/EM)",
"CD4+ Naive T cells", "gamma-delta T cells",
"Late prolif: CCNB1/2+ CDK1+ T cells",
"Early prolif: MCM3/5/7+ PCNA+ T cells",
"Late prolif: histones enriched MKI67+ T cells",
"Late prolif: STMN1+ BIRC5+",
"Ribosomal/Mitochondrial/Degraded cells"))
occRep$cloneSize <- factor(occRep$cloneSize, levels = c("Small (1 < X <= 5)", "Medium (5 < X <= 20)", "Large (20 < X <= 100)", "Hyperexpanded (100 < X <= 500)"))
# Draw barplot with grouping & stacking
ggplot(occRep, aes(x = Timepoint, y = n, fill = cloneSize)) +
geom_bar(stat = "identity", position = "stack") +
facet_grid(~ annotation) +
scale_fill_manual(values = col_clono) + facet_grid2(~ annotation, strip = strip) +
theme(strip.text = element_text(colour = NA),
axis.text.x = element_text(angle = 45, vjust = 0.5, hjust=0.5))
Plot Fig3S D
I <- getCirclize(data_pt1_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 1")
P <- getCirclize(data_pt1_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 1")
I <- getCirclize(data_pt2_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 2")
P <- getCirclize(data_pt2_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 2")
I <- getCirclize(data_pt3_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 3")
P <- getCirclize(data_pt3_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 3")
I <- getCirclize(data_pt4_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 4")
P <- getCirclize(data_pt4_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 4")
I <- getCirclize(data_pt5_I, group.by = "annotation")
I<-chordDiagram(I, self.link = 1, grid.col =cols_anno , annotationTrack = c("grid", "axis"))
title("IP Patient 5")
P <- getCirclize(data_pt5_P, group.by = "annotation")
P<-chordDiagram(P, self.link = 1, grid.col =cols_anno, annotationTrack = c("grid", "axis") )
title("Peak Patient 5")
Plot Fig3S E
text <- element_text(color = "black", size = 16)
text.lab <- element_text(color = "black", size = 12)
text.lab2 <- element_text(color = "black", size = 10)
text.lab3 <- element_text(color = "black", size = 9)
aux_df <- data.frame(barcodes=colnames(data2),idents=data2$annotation,Pseudotime=data2$Pseudotime, Timepoint = data2$Timepoint)
ident_trajectory <- ggplot(aux_df,aes(x=idents,y=Pseudotime,col=Timepoint)) +
geom_quasirandom(alpha=.8) + coord_flip() + theme_classic() +
scale_color_manual(values=cols_timepoint) + theme(text = element_text(size = 14))
ident_trajectory
Plot Fig3S F
ordering.genes <- VariableFeatures(seurat)
marker_genes <- row.names(subset(fData(data2),
gene_short_name %in% ordering.genes))
data2 <- reduceDimension(data2,
max_components = 2,
norm_method = 'log',
num_dim = 3,
reduction_method = 'tSNE',
verbose = T)
data2 <- clusterCells(data2, verbose = F, num_clusters = 10)diff_test_res <- differentialGeneTest(data2[marker_genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime)",
cores=4)
sig_gene_names <- row.names(subset(diff_test_res, qval < 0.000000000000000000000000000000000000001))sig_gene_names <- unique(c("STMN1", "CDC20", "HIST2H2BF", "UBE2C", "GNLY", "GZMK", "CCNB1", "MKI67", "KLRG1", "KLRB1", "CD69", "GZMH", "HMGB1", "PCNA","GZMB", "IFITM1", "CD27", "LAG3", "MKI67", "TOP2A", "CDK1", "CXCR3", "GZMK", "NKG7", "HSPA1A", "HSPA1B", "PRF1", "TRDV1", "TRGV2" ))
plot_pseudotime_heatmap(data2[sig_gene_names,],
num_clusters = 3,
cores = 1,
show_rownames = T,
hmcols = colorRampPalette(colors = c("blue", "white", "red"))(62))
Plot Fig3S G
cols2=c("CD4+ CEntral/Effector memory T cells (CM/EM)"="white",
"CD4+ Naive T cells"="white",
"CD8+ cytotoxic T cells"="#185B88",
"CD8+ Eff/Mem T cells (EM)"="white",
"CD8+ Effector T cells (E)"="white",
"Early prolif: MCM3/5/7+ PCNA+ T cells"="white",
"Late prolif: CCNB1/2+ CDK1+ T cells"="white",
"Late prolif: histones enriched MKI67+ T cells"="white",
"Late prolif: STMN1+ BIRC5+"="white",
"Ribosomal/Mitochondrial/Degraded cells"="white")
data_pt1 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient1")
data_pt2 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient2")
data_pt3 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient3")
data_pt4 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient4")
data_pt5 <- subset(seurat_carneg_nogd, subset = Patient_id == "patient5")
pt1 <-alluvialClones(data_pt1, cloneCall = "CTstrict",
y.axes = c("annotation", "cloneType", "Timepoint"),
color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 1650))
pt2 <-alluvialClones(data_pt2, cloneCall = "CTstrict",
y.axes = c("annotation", "cloneType", "Timepoint"),
color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 1650))
pt3 <-alluvialClones(data_pt3, cloneCall = "CTstrict",
y.axes = c("annotation", "cloneType", "Timepoint"),
color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 1650))
pt4 <-alluvialClones(data_pt4, cloneCall = "CTstrict",
y.axes = c("annotation", "cloneType", "Timepoint"),
color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") +
scale_y_continuous(limits = c(0, 1650))
pt5 <-alluvialClones(data_pt5, cloneCall = "CTstrict",
y.axes = c("annotation", "cloneType", "Timepoint"),
color = "annotation") + scale_fill_manual(values=cols2)+theme(legend.position = "none") + scale_y_continuous(limits = c(0, 1650))
ggarrange(pt1 + facet_grid(. ~"Patient 1"),
pt2 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 2"),
pt3 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 3"),
pt4 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 4"),
pt5 +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
axis.title.y = element_blank(),
plot.margin = margin(r = 1, l = 1) )+ facet_grid(. ~"Patient 5"), align = "hv", ncol=5)
Plot Fig3S H
seurat_carpos <- subset(seurat, subset = Class1 == "CAR+")
table_unique_clono <- scRepertoire:::.expression2List(seurat_carpos, split.by = "Sample_id")
table_unique_clono <-clonalQuant(table_unique_clono, exportTable = T)
table_unique_clono$patient <- revalue(table_unique_clono$values, c("patient1_IP"="Patient 1", "patient1_Peak"="Patient 1", "patient2_IP"="Patient 2", "patient2_Peak"="Patient 2", "patient3_IP"="Patient 3", "patient3_Peak"="Patient 3", "patient4_IP"="Patient 4", "patient4_Peak"="Patient 4", "patient5_IP"="Patient 5", "patient5_Peak"="Patient 5" ))
table_unique_clono$patient <- factor(table_unique_clono$patient)
table_unique_clono$timepoint <- rep(c("IP","Peak"),5)
ggplot(table_unique_clono, aes(fill=values, y=contigs, x=patient)) +
geom_bar(position="dodge", stat="identity")+ scale_fill_manual(values=cols_patient2) +
geom_text(aes(label=timepoint), position = position_dodge(width = .9), vjust = -0.3, color = "black") + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))
Session Info
sessionInfo()## R version 4.2.0 (2022-04-22)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Big Sur/Monterey 10.16
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] splines stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] monocle_2.24.0 DDRTree_0.1.5 irlba_2.3.5.1 VGAM_1.1-9 Biobase_2.58.0 BiocGenerics_0.44.0 Matrix_1.6-5 ggbeeswarm_0.7.2 ggh4x_0.2.8 dplyr_1.1.4 plyr_1.8.9 stringr_1.5.1 Startrac_0.1.0 circlize_0.4.16 scRepertoire_2.0.4 cowplot_1.1.3 ggpubr_0.6.0 ggplot2_3.5.1 reshape2_1.4.4 Seurat_5.0.1 SeuratObject_5.0.1 sp_2.1-2 gdata_3.0.0
##
## loaded via a namespace (and not attached):
## [1] utf8_1.2.4 spatstat.explore_3.2-5 reticulate_1.34.0 tidyselect_1.2.1 htmlwidgets_1.6.4 docopt_0.7.1 combinat_0.0-8 grid_4.2.0 Rtsne_0.17 munsell_0.5.1 codetools_0.2-20 ica_1.0-3 future_1.33.2 miniUI_0.1.1.1 withr_3.0.0 fastICA_1.2-4 spatstat.random_3.2-2 colorspace_2.1-0 progressr_0.14.0 highr_0.10 knitr_1.41 ggalluvial_0.12.5 rstudioapi_0.16.0 SingleCellExperiment_1.20.1 ROCR_1.0-11 ggsignif_0.6.4 tensor_1.5 listenv_0.9.1 labeling_0.4.3 MatrixGenerics_1.10.0 slam_0.1-50 GenomeInfoDbData_1.2.9 polyclip_1.10-6 pheatmap_1.0.12 farver_2.1.1 parallelly_1.36.0 vctrs_0.6.5 generics_0.1.3 xfun_0.41 doParallel_1.0.17 R6_2.5.1 GenomeInfoDb_1.34.9
## [43] clue_0.3-65 graphlayouts_1.0.1 bitops_1.0-7 spatstat.utils_3.0-4 cachem_1.0.8 DelayedArray_0.24.0 promises_1.2.1 scales_1.3.0 ggraph_2.1.0 beeswarm_0.4.0 gtable_0.3.5 globals_0.16.3 goftest_1.2-3 spam_2.10-0 tidygraph_1.3.0 MatrixModels_0.5-3 rlang_1.1.2 GlobalOptions_0.1.2 rstatix_0.7.2 lazyeval_0.2.2 spatstat.geom_3.2-7 broom_1.0.5 yaml_2.3.8 abind_1.4-5 backports_1.4.1 httpuv_1.6.13 tools_4.2.0 cubature_2.1.0 jquerylib_0.1.4 RColorBrewer_1.1-3 ggdendro_0.2.0 ggridges_0.5.6 hash_2.2.6.3 Rcpp_1.0.11 zlibbioc_1.44.0 purrr_1.0.2 RCurl_1.98-1.13 deldir_2.0-2 GetoptLong_1.0.5 pbapply_1.7-2 viridis_0.6.5 S4Vectors_0.36.2
## [85] zoo_1.8-12 iNEXT_3.0.1 SummarizedExperiment_1.28.0 ggrepel_0.9.4 cluster_2.1.6 magrittr_2.0.3 data.table_1.15.4 RSpectra_0.16-1 scattermore_1.2 SparseM_1.81 lmtest_0.9-40 RANN_2.6.1 truncdist_1.0-2 fitdistrplus_1.1-11 matrixStats_1.2.0 gsl_2.1-8 patchwork_1.2.0.9000 mime_0.12 evaluate_0.23 xtable_1.8-4 sparsesvd_0.2-2 shape_1.4.6.1 fastDummies_1.7.3 IRanges_2.32.0 gridExtra_2.3 HSMMSingleCell_1.18.0 compiler_4.2.0 tibble_3.2.1 crayon_1.5.2 KernSmooth_2.23-22 htmltools_0.5.7 later_1.3.2 tidyr_1.3.0 tweenr_2.0.2 ComplexHeatmap_2.14.0 MASS_7.3-57 leidenbase_0.1.27 car_3.1-2 cli_3.6.2 evd_2.3-6.1 parallel_4.2.0 dotCall64_1.1-1
## [127] igraph_1.5.1 GenomicRanges_1.50.2 pkgconfig_2.0.3 plotly_4.10.4 spatstat.sparse_3.0-3 foreach_1.5.2 vipor_0.4.7 bslib_0.4.1 stringdist_0.9.12 XVector_0.38.0 digest_0.6.33 sctransform_0.4.1 RcppAnnoy_0.0.21 spatstat.data_3.0-4 rmarkdown_2.18 leiden_0.4.3.1 uwot_0.1.16 evmix_2.12 shiny_1.8.1.1 gtools_3.9.5 quantreg_5.97 rjson_0.2.21 lifecycle_1.0.4 nlme_3.1-164 jsonlite_1.8.8 carData_3.0-5 limma_3.54.2 viridisLite_0.4.2 fansi_1.0.6 pillar_1.9.0 lattice_0.22-5 fastmap_1.1.1 httr_1.4.7 survival_3.5-7 glue_1.6.2 qlcMatrix_0.9.7 iterators_1.0.14 png_0.1-8 ggforce_0.4.1 stringi_1.8.3 sass_0.4.8 RcppHNSW_0.5.0
## [169] future.apply_1.11.2